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ABSTRACT 


Viscous flow past a circular cylinder becomes unstable around 
Reynolds number Re = 40. With a numerical technique based on Newton's 
method and made possible by the use of a supercomputer, steady (but 
unstable) solutions have been calculated up to Re = 400. It is found 
that the wake continues to grow in length approximately linearly with 
Re. However, in conflict with available asymptotic predictions, the 
width starts to increase very rapidly around Re = 500. All numerical 
calculations have been performed on the CDC Cyber 205 at the CDC 
Service Center in Arden Hills, Minnesota. 


INTRODUCTION 

The structure of viscous steady flow past a circular cylinder at 
high Reynolds numbers forms one of the classical problems in fluid 
mechanics. In spite of much attention, several fundamental questions 
remain open. Apart from a previous calculation by the present author 
[6], complete, steady flow fields have been obtained numerically only 
up to around Re = 100. This is also close to the upper limit for 
experiments (due to temporal instabilities). Both the early numerics 
and the experiments point to a recirculation region growing linearly 
in length with Re. Figure 1 shows the length of the wake bubble 
against Reynolds number according to some different calculations. 
Persistence of this growth for Re -> oo has been assumed in most 
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recent asymptotic studies of steady high Reynolds number flows past a 
body (e.g. F.T. Smith [ 1 3 ] ) - A possible Euler flow, consistent with 
this idea, was analyzed by Brodetsky [3] in 1923. It is known as the 
Helmholtz- Kirchhoff free streamline model. This suggested limit is 
characterized by two vortex sheets leaving the body tangentially 
approximately 55° from the upwind center line and extending to 
downstream infinity, enclosing a region of stagnant flow. Although 
this undoubtedly is a solution for Re = 00, G.K. Batchelor [2] gave 
in 1956 several arguments against this being a possible limit for 
Re -> 00. He proposed an alternative in which a finite wake with 
piecewise constant vorticity was bounded by vortex sheets. Some 
suggestions about how such a flow might be reached as a limit for 
increasing Reynolds number have been given by Peregrine [lO], 
However, only very few Euler solutions of this so called 
Prandtl-Batchelor type have been calculated (e.g.[l2] contains one 
example and some further references). None of these are for flow past 
a cylinder. Figure 2 gives an 'artists impression' of what the two 
models for infinite Re might look like. The calculation [6] hinted 
at a process leading to a shortening of the wake. The present work 
suggests (in agreement with F.T. Smith [ 1 4 ] ) this shortening at 
Re = 300 was erroneous and caused by insufficient numerical 
resolution. However, our best current evidence is that the 
qualitative result was correct. We beleive that a reversal of trends 
towards a shorter wake can be expected around Re = 500. This 
contrasts with the conclusions in [14]. Our main evidence is that 
the wake increases in width far more rapidly after Re = 300 than the 
asymtotic analysis allows for. Independently of the position of 
artificial boundaries and of numerical resolution, we find that the 
flow is of different character past Re = 300. Significant amounts of 
vorticity are then re-circulated back into the wake bubble from its 
end. We hope to soon carry this study past Re = 400. 

All the numerical calculations in this present work were 
performed on the Control Data Corporation Cyber 205 computer located 
at the CDC Service Center in Arden Hills, Minnesota. We wish to 
express our gratitude to Control Data Corporation for making this 
system available for this work. 
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MATHEMATICAL FORMULATION. 


With a cylinder of radius 1 and a Reynolds number based on the 
diameter, the governing time independent Navier-Stokes equations, 
expressed in streamfunction and vorticity oq , take the form: 
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Accurate numerical approximation and economical computational 
solution of these equations in the given geometry poses a series of 
difficulties which previous investigators have dealt with in a 
variety of ways. The most serious of the difficulties seem to be: 


1 . Boundary comditions for at large distances. 

2. Boundary condition for U> at the body surface. 

Avoiding the loss of accuracy that comes with upwind 
differencing. 

4. Economical choice of computational grid. 

5* Reliable and fast rate of convergence of numerical 
iterations. 

The point 5 above has been the limiting factor in virtually all 
previous attempts to reach high Reynolds numbers. No reliable 
technique has emerged to prevent slowly converging iteration schemes 
from picking up physical instabilities in the artificial time of the 
iterations . 
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NUMERICAL METHOD 


All vorticity is concentrated on the body surface and in a quite 
thin wake downstream of the body. Outside this region we can use the 
much simpler equations: 


(3) 

(4) 


= o 


CO 


= 0 


The top part of Figure 3 shows the upper half plane minus a unit 
circle and, dotted, a region which contains all the vorticity (apart 
from the far wake). The bottom part of the figure shows how the 
mapping z = t /*•+ i/iyr maps these to the first quadrant and a 
rectangle respectively. Figure 4 shows what a rectangular grid in 
the z-plane (with non-uniform stretching in the vertical direction) 
can look like in the x-plane. The Navier-Stokes equations, 


transformed to 
(2): 

the z- 

■plane 

take a 

form almost identical to (l) and 
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modified further 

is the Jacobian 
by subtracting 

of the mapping. These equations 
out potential flow. The stream 


function for the difference is r - x -2 . On a grid in the 

(stretched) z-plane, equations (5) and (6) were approximated at all 
interior points with centered second order finite differences. To 
close the system, boundary conditions have to be implemented for 
and to at all boundaries. 


The extreme sensitivity of the final solution to small errors in 
these conditions has only recently been fully recognized [6]. For 
example already at Re = 2 it was found that use of the free stram 
value for along circular outer boundaries at distances 23.1 and 
91.5 caused 18 % and 4.4 % errors in the level of vorticity on the 
body surface. 
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The ’Oseen' approximation is the leading term in an asymptotic 
expansion for the flow far out in a wake (e.g. Imai [8]). In polar 
coordinates, it takes the form 

C 6 

(7) 4^ = ( “ erf Q ) 


C^Re Q -Q 

(8) to = e 

4 r 

where Q = (^Re r) /Z sin-=^0 , erf Q = 2’K" ^ ds and C & the drag 

coefficient. C p can be evaluated as a line integral around the body. 

The performance of this Oseen condition as an outer boundary 
condition is disappointing. The percentage errors mentioned above 
improve, but only to to 3-4 % and 1.2 % respectively. For increasing 
Re, direct use of (7) becomes meaningless. Figure 5 illustrates this 
by comparing the true 4* (here the difference between streamfunction 
and free stream, not potential flow) with the values from (7) at 
Re = 200. The two fields bear no resemblance to each other at the 
distances from the body we are interested in. 

Comparison with numerics suggest that (8) is far more accurate 
than (7). Furthermore 

1 . Any errors in (8) are present only in a very narrow region 
along the outflow axis, not along the whole upper boundary 
as with (7). 

2. The governing equation for lO is of a type which cannot 
transport incorrect information for Ca 3 back up towards the 
cylinder. 

With this background, let us briefly outline how the boundary 
conditions of high accuracy can be implemented on the edges of the 
present computational region. Figure 6 shows this region in the 
z-plane with a typical vorticity field together with its reflections 
in the coordinate axis. 
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BOUNDARY CONDITIONS FOR 4\ 

Left boundary: = 0 ,0 <f)< 9vj* 

Bottom boundary: 0 ,0 < 

Right boundary: J=S M ,0 


Top boundary: *)- 0,^0 <"£<5^. 

A correction to the integral above for vorticity reaching 
outside the downstream boundary can easily be incorporated. For a 
fixed grid, the dependence of 4" at each boundary point on io at 
each internal point is independent of Re and can be calculated as a 
large matrix once and for all. A boundary condition of this kind was 
used in all the calculations presented below. However, we currently 
use a different condition. A wide two level difference formula can be 
found which is consistent only with the decaying modes of the 
equation £ 4* - 0 (as opposed to the usual 5-point 5-level formula 
used inside the region to approximate both growing and decaying 
modes) . 


= 0. 


4>= 0. 

Jj\v 

- to (noting that << 

<< ^“ 2 . along this boundary). 
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BOUNDARY CONDITIONS FOR U3 . 


Left boundary: 0 ,0 <*)<>()*> . 

Bottom boundary: 0 ,0 < £< 2 . 

2 i s<t M . 

Right boundary: l=K’° 

Top boundary: ^=9*1,0 ^ ? <S M . 


<0 = 0 . 

A relation based on 0 

and an even function of n . 

10 = 0 . 

Ok (5/s ) 2 

to = o. 


The condition at the right boundary comes from the observation 
that the leading term of (8), transformed to 5 -coordinates 
simplifies to 


(9) W 



-c*9 


where c , and c x are constants. The mapping has achieved a 
separation of variables. 
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The discrete approximations at the interior points together with 
the boundary conditions form, after minor simplifications (explicitly 
eliminating all boundary unknowns apart from ^ at the top 
boundary), a non-linear algebraic system of (M-2)(2N-3) equations 
with equally many unknowns. In most earlier works, great care has 
been taken to ensure that, at this stage, this (or some equivalent) 
non-linear system has a diagonally dominant form for low Re. This 
would allow direct functional iteration to convergence. Techniques 
like upwind differencing [l],[4],[ll] help in this respect at the 
cost of lowered accuracy. Newton's method, described below, offers an 
outstanding alternative. 


NEWTON'S METHOD. 


Newton's method is a very well known procedure for finding zeros 
of scalar functions. If a function f(x) is given, we can find an x 
such that f(x)-0 by the procedure: 

(10) x 0 'close' guess of root 

f(x J 

n - 0, 1 , 2, ... 

f(x n ) 

The iteration step can be written 

(12) f'(x^) -A x ^ = ”f(x H ) 

Known, f' evalu- 
ated at the latest 
available approxi- 
mation x ^ . 


Written in this form, the generalization to systems is 
straightforward. For example the system with three equations in 
three unknowns: 


Unknown, the 
correction we 
should apply 
to x n , i. e. 

%,r V Ax n- 


Known, residual. 
Should be zero 
if x n had been 
exact. 


( 11 ) 
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( 15 ) 


f(x, y, z) = 0 

g(x, y, z) = 0 

h(x, y, z) = 0 


can be iterated 
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LA* j 

L h(x ,y ,z 
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Known, 

"Jacobian" 

Unknown, 

Known, 



of system. 


corrections. 

residual. 



Each iteration 

involves 

the 

solution of a linear system. 

Like in 

the 

scalar case, convergence is 

quadratic and 

guaranteed 

to occur 

for 

approximations 

sufficiently 

close to any 

simple solution. 

The 


realization that this procedure is practical for extremely large 
systems (several thousands of equations) is rather recent and linked 
to the emergence of powerful computers. 

For our present problem, use of Newton's method offers several 
major advantages: 

1. The quadratic convergence allows no possibility of 'inheriting' 
temporal instabilities to the artificial time of the iterations. 
Convergence is guaranteed if an isolated solution exists in the 
neighborhood of a guess. 

2. If turning points or bifurcation points are found, they will 
cause no difficulties. 

3* No upwind differencing is needed. This procedure is typically 
employed for two reasons: 
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1 . To ensure convergence of an iterative method. 

2. To avoid mesh size oscillations. 

The first reason no longer applies. The second one alone can 
then be addressed in more refined ways. 

4. Boundary conditions at the body surface become easier to 
implement. The fact that we have two conditions on W and none 
on c-o can cause a problem if (5) and (6) are treated separately. 
With Newton's method, all we need is that the number of 
conditions is right. 

The only disadvantage with Newton's method is the computational 
cost. This is where supercomputers enters our picture. 


SOLUTION 0? LINEAR SYSTEM 

Let [^J., j=2,3,...,N be vectors with M* -values from grid lines 
2,3,... ,11 and similarly for <0^ (j=2,..., N-l). For example would 
contain the -values along the grid row nearest to the 5 -axis and 
the values along the top boundary. The structure of the entries in 
the Jacobian matrix reflects directly on the difference stencils and 
the boundary conditions. Figure 7 shows a suitable ordering of 
equations and unknowns and the corresponding structure of the 
Jacobian. Since the top right corner contains a single diagonal, 
explicit multiples of the top (N-2)(M-2) equations can be superposed 

on the equations below to modify the structure to the one in Figure 

8. The bottom left corner form a separated system of size (N-l)(M-2). 
This system was solved by a border algorithm similar to the one 
described in [ 9 ]. The major cost comes from the LU-factorization of 
A. However, one more rearrangement can be done to achieve a 
significant speedup. The A-matrix has a block 5-diagonal form with 
the structure shown in Figure 8. A similarity transform with a 
permutation matrix can rearrange this into another matrix of 
identical structure. Instead of N-2 rows of blocks, each of size M-2 , 
we get M-2 rows of blocks of size N-2. With M typically around 6*N 
and cost proportional to the square of the bandwidth, this reduces 

the memory needed for the LU-decomposition about 6 and the operation 

count by 36. 

The complete linear solver lends itself ideally to 
vectorization. Every part of significant cost turns out take a form 
of a 'linked triad* with vectors never shorter than 4(N-2)+1 or M. 
The linked triad on the Cyber 205 is the fastest floating point 
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operation the machine offers. Expressions of the form 
vector-op-vector-op-scalar where one ’op' is + or -, the other * can 
execute with both operations running simultaneously. On the 2-pipe 
205, the algorithm has a potential for 200 mflops (million floating 
point operations per second, 64-bit accuracy). Including a startup 
cost of 83 macihne cycles per linked triad operation, average vector 
length of around 166 (which we will exceed in later test cases) could 
give a full 100 mflop overall computational performance. In the 
calculations presented below, the grid had 131 by 21 points. 
Building up the Jacobian (in scalar mode) takes 2.3 seconds and the 
solution of the linear system 3*7 seconds (for an average of 55 
mflops during this part). Recently implemented vectorization of the 
Jacobian and the new boundary condition brings these numbers to 0.026 
seconds, 1.75 seconds and 60 mflops respectively. 


PHYSICAL CONCLUSIONS 


This report is a preliminary one of work in progress. Only a few 
initial test runs have been performed so far. However, we can already 
conclude that the wake appears to continue a linear growth in length 
with increasing Reynolds numbers up to Re = 400. Figure 9 shows wake 
length versus Re for some previous calculations compared with current 
results. Figure 10 shows streamlines and Figure 11 vorticity fields 
for different values of Re up to 4 00. The vorticity field at Re = 
400 shows a recirculation back into the wake from the end of the 
bubble as well as a quite sudden increase in width. Our most recent 
tests with a computational grid of 196 by 31 points (density 
increased by 3/2 in each direction) leaves these features completely 
unchanged. The onset occurs near Re = 300 and the widening progresses 
at a rate which can be determined accurately and which far exceeds 
the one predicted by available asymptotic models. 

The flow fields in figures 10 and 11 were obtained from a 131*21 
grid in the z-plane with 
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(14) 

(15) 


, i=0 , 1 , ... ,150 


5; = i/12 

3 

n. =«.£.+ (1-4*)C, Vr- j/18 , J-0,1, ... ,20 , K-0.15 

) 3 4 4 

This places the right boundary at a distance 115*4 from the center of 
the cylinder. Preliminary tests involving moving this and the top 
boundaries in and out suggest that they are sufficiently far out with 
the present choice of grid. Figure 4 showed part of this grid. 

The major open questions at the moment are: 

Physically: 

1 . Will the wake keep on growing? 

2. Are there any other branches of solutions (bifurcations 
etc . ) ? 

Numerically: 

1. Is there any alternative to Newton's method which still 
possesses a reliable rate of convergence? 

2. Is there any faster way than Gaussian elimination to solve 
the linear system in Newton's method? 

At present, the numerical questions are wide open and of 
fundamental importance to many other applications as well. Current 
numerical methods together with vector computers like the Cyber 205 
probably form sufficiently powerful tools to settle conclusively the 
physical questions raised here. 
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Fornberg (1980) 
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Figure L. Length of wake bubble for low Reynolds numbers according to some different calculations 



Figure 2. Schematic illustration of free streamline and 
the Prandtl-Batchelor models . 
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Figure 6. Typical vorLicity field in complex z-plane 
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Line with slope -^=0.17^ 



Figure 9: Length of wake bubble for different Reynolds numbers. 
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Reynolds numbers. 




Contours of constant vorticity (values w = -12., 
0., .1, .2, .5) for different Reynolds numbers. 
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